rm(list=ls())
set.seed(8768)
library(fGarch)
library
library(MASS)
#### Output file
nom="Backtest-VaR-exchgrate-insample.out"

##### Load the file */
garch = read.csv("exchg.csv",header = FALSE, sep =";")
nmax=ncol(garch)
n=nrow(garch)-1

### log p.d.f.
vrstu<-function(x,nu){ 
if (nu<0){likeli=rep(-1e8,length(x))}
else{lknu=log(gamma((nu+1)/2))-0.5*log(pi*(nu-2))-log(gamma(nu/2));
likeli=lknu-((nu+1)/2)*log(1+(x^2)/(nu-2));
}
return(likeli)}


### Score Standardized student

scorestu<-function(x,nu){
knu=0.5*(-1/(nu-2)+digamma((nu+1)/2)-digamma(nu/2));
fst=(nu+1)/(2*(nu-2))*((x^2)/(nu-2+x^2))-log(1+(x^2)/(nu-2))/2
return(fst+knu)}



    ##Log-likelihood
    garchLLH = function(x,parm) {
        omega = parm[1]; alpha = parm[2]; beta = parm[3];nu=parm[4];
        z =x; ltv = mean(z^2)
        # Use Filter Representation:
        e = omega + alpha * c(ltv, z[-length(x)]^2)
        h = filter(e, beta, "r", init =ltv )
        llh = sum(0.5*log(h)-vrstu(z/sqrt(h),nu))/length(x);
        return(llh)
    }

    #### Log likekihod with constraint  on the parameters
    vrais=function(x,b){
		bb=b;
		bb[1]=b[1]^2;
		bb[3]=exp(-b[2]^2-b[3]^2);
		bb[2]=exp(-b[2]^2)-bb[3];
		return(garchLLH(x,bb))
	}
stat_exchg=function(num,pp){
u0=garch[,num]
data=log(u0[-1]/u0[1:n])####log-returns
	#Initial conditions
    	data_est=data
	vol=var(data_est)
    	adeb=0.1
    	bdeb=0.8
    	odeb=vol*(1-adeb-bdeb);
	tdeb=c(sqrt(odeb),sqrt(-log(adeb+bdeb)),sqrt(log((adeb+bdeb)/bdeb)),8)

#######MLE 
    
	log_lik=function(param){vrais(data_est,param)}

    options=list("algorithm"="NLOPT_LN_NELDERMEAD","xtol_rel"=1e-8,"maxeval"=1000000)
	res<-nloptr(x0=tdeb,eval_f=log_lik,opts=options)
	estimpar=res$solution


	omega=estimpar[1]^2
	beta=exp(-estimpar[2]^2-estimpar[3]^2);
	alpha=exp(-estimpar[2]^2)-beta  
    	nu=(estimpar[4])
   valest=c(omega,alpha,beta,nu)


    z =data; 
    ltv = mean(z^2)
    e = omega + alpha * c(ltv, z[-length(z)]^2)
    h = filter(e, beta, "r", init =ltv )
    c3=filter(c(0,h[-n]),beta,"r",init=0)
    c2=filter(c(0,z[-n]^2),beta,"r",init=0)
    h=as.vector(h)
    ### The derivative of log vol with respect to the parameters.
    der=cbind(1/h/(1-beta),c2/h,c3/h)
    der[1,]=c(0,0,0)

   ######### Hit sequence for alpha=pp
    	epst=as.vector(as.numeric((data/sqrt(h))))
	#### Hits
    	It=1*(pt(epst*sqrt(nu/(nu-2)),nu)<pp)
    	qa=qt(pp,nu)/sqrt(nu/(nu-2))
	coef=qt(pp,nu)*dt(qt(pp,nu),nu)
      ### Correction wrt the true score function
	sc=((nu+1)*(epst^2)/(nu-2+epst^2)-1)
    	fulsc=cbind(sc,scorestu(epst,nu))
    	matV=(t(fulsc)%*%fulsc/n)
    	matV[1,1]=(2*nu)/(nu+3)
    	matP=as.matrix(apply((It-pp)*fulsc,2,mean))
    	matP[1]=-coef #This is MINUS Expectation of "d It/dtheta"
    	dlog=der/2
    	espder=apply(dlog,2,mean)
    	Ahat=cbind(matP[1]*t(espder),matP[2])
    	score=cbind(dlog*sc,scorestu(epst,nu))
    	Vhat=ginv(t(score)%*%score/n)
    	corscore=(Ahat)%*%Vhat%*%t(Ahat)
	############### Corrected hits ###########################
	###### e_t in the paper, i.e. projection onto the orthogonal of an auxiliaray score
    	et=It-pp-fulsc%*%ginv(matV)%*%matP
    	# Normalization to have unit variance
    	et=et/as.numeric(sqrt(pp*(1-pp)-t(matP)%*%ginv(matV)%*%matP))

    	##### etstar in the paper, i.e. projection onto the orthogonal of the true score
    	etstar=It-pp-score%*%Vhat%*%t(Ahat)
    	# Normalization to have unit variance
    	etst=etstar/as.numeric(sqrt((pp*(1-pp)-corscore)))

    	##### Correction of the variance of the hit sequence
    	Itbis=(It-pp)/as.numeric(sqrt((pp*(1-pp)-corscore)))

	
	##### Unconditional tests
      xiM=n*mean(et)^2
     	# Project orthogonal score 
    	xiMb=n*mean(etst)^2
    	#Correction paruncertainty 
    	xiPU=n*mean(Itbis)^2
    	#Hit without correcting
    	xihit=n*mean(It-pp)^2/as.numeric(pp*(1-pp))
    	
	#### Joint moments e_t.e_{t-h}
    	etm1=c(NA,et[-length(et)])
    	etm2=c(NA,etm1[-length(et)])
    	etm3=c(NA,etm2[-length(et)])
      xiet1=(n-1)*mean(etm1*et,na.rm=TRUE)^2
    	xiet2=(n-2)*mean(etm2*et,na.rm=TRUE)^2
    	xiet3=(n-3)*mean(etm3*et,na.rm=TRUE)^2

	#### Joint moments est_t.est_{t-h}
    	etstm1=c(NA,etst[-length(etst)])
    	etstm2=c(NA,etstm1[-length(etst)])
    	etstm3=c(NA,etstm2[-length(etst)])
      xietst1=(n-1)*mean(etstm1*etst,na.rm=TRUE)^2
    	xietst2=(n-2)*mean(etstm2*etst,na.rm=TRUE)^2
    	xietst3=(n-3)*mean(etstm3*etst,na.rm=TRUE)^2

	##### Joint moments I_t.I_{t-h} corrected for par. uncertainty
    	Itloc=It-pp	
    	Itm1=c(NA,Itloc[-length(It)])
    	Itm2=c(NA,Itm1[-length(It)])
    	Itm3=c(NA,Itm2[-length(It)])
      xihit1=(n-1)*mean(Itm1*Itloc,na.rm=TRUE)^2
    	Ahat1=c(mean(Itloc*Itm1*score[,1],na.rm=TRUE),mean(Itloc*Itm1*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm1*score[,3],na.rm=TRUE),mean(Itloc*Itm1*score[,4],na.rm=TRUE))
    	corscore1=t(Ahat1)%*%Vhat%*%Ahat1
    	xihit1c=xihit1/(pp*(1-pp)*pp*(1-pp)-corscore1)
    	xihit1sc=xihit1/(pp*(1-pp)*pp*(1-pp))
    	

    	xihit2=(n-2)*mean(Itm2*Itloc,na.rm=TRUE)^2
    	Ahat2=c(mean(Itloc*Itm2*score[,1],na.rm=TRUE),mean(Itloc*Itm2*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm2*score[,3],na.rm=TRUE),mean(Itloc*Itm2*score[,4],na.rm=TRUE))
    	corscore2=t(Ahat2)%*%Vhat%*%Ahat2
    	xihit2c=xihit2/(pp*(1-pp)*pp*(1-pp)-corscore2)
    	xihit2sc=xihit2/(pp*(1-pp)*pp*(1-pp))
    	
	xihit3=(n-3)*mean(Itm3*Itloc,na.rm=TRUE)^2
    	Ahat3=c(mean(Itloc*Itm3*score[,1],na.rm=TRUE),mean(Itloc*Itm3*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm3*score[,3],na.rm=TRUE),mean(Itloc*Itm3*score[,4],na.rm=TRUE))
    	corscore3=t(Ahat3)%*%Vhat%*%Ahat3
    	xihit3c=xihit3/(pp*(1-pp)*pp*(1-pp)-corscore3)
    	xihit3sc=xihit3/(pp*(1-pp)*pp*(1-pp))
    	

	testavecnu=c(xiM,xiMb,xiPU,xihit,xiet1,xiet2,xiet3,xietst1,xietst2,xietst3,xihit1c,xihit2c,xihit3c,xihit1sc,xihit2sc,xihit3sc)


	return(c(valest,(1-pchisq(c(testavecnu),1)),mean(It)))
}


######OUTPUT#####################

outputres<-matrix(rep(NA,4*nmax),ncol=nmax)
outputhit<-matrix(rep(NA,4*nmax),ncol=nmax)
outputtableavecnu<-matrix(rep(NA,16*(3*nmax+3)),ncol=nmax*3+3)

colnames(outputtableavecnu) <-c(" ","UK-US","FF-US","SF-US","Yen-US"," ","UK-US","FF-US","SF-US","Yen-US"," ","UK-US","FF-US","SF-US","Yen-US")

rownames(outputtableavecnu) <-c("$e_t$","$e_t^*$","$e_t^c$","$I_t$",
"$e_te_{t-1}$","$e_te_{t-2}$","$e_te_{t-3}$",
"$e_t^*e_{t-1}^*$","$e_t^*e_{t-2}^*$","$e_t^*e_{t-3}^*$","$e_t^ce_{t-1}^c$","$e_t^ce_{t-2}^c$","$e_t^ce_{t-3}^c$",
"$I_tI_{t-1}$","$I_tI_{t-2}$","$I_tI_{t-3}$")

listep=c(0.005,0.01,0.05)
ss=1
is=1
for(pp in listep ){
	ss=ss+1
	is=is+1
	for (numcol in 1:nmax){
	output=stat_exchg(numcol,pp)
	outputres[,numcol]=output[1:4]
	outputtableavecnu[,ss]=output[5:20]
	outputhit[is,numcol]=output[21]
	ss=ss+1
}
}
library(xtable)
outputfile=paste(nom,".tex")
print(xtable(outputtableavecnu,digits=3),sanitize.text.function=function(x){x},hline.after = getOption("xtable.hline.after", c(3,6,9)),file=outputfile)
print(xtable(outputres,digits=3))
print(xtable(outputhit,digits=3))
